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We study, using Monte Carlo simulations, the interaction between infinite heterogeneously charged surfaces 
inside an electrolyte solution. The surfaces are overall neutral with quenched charged domains. An average 
over the quenched disorder is performed to obtain the net force. We find that the interaction between the 
surfaces is repulsive at short distances and is attractive for larger separations. 


I. INTRODUCTION 

In physical chemistry and biophysics one often finds 
situations in which electrolyte solution is confined be¬ 
tween charged surfaces. The surfaces can belong to 
macromolecules, colloidal particles, electrodes or mem¬ 
branes. Presence of electrolyte between surfaces can 
strongly modify the interaction between themjii^. 

Over seventy years ago, Derjaguin, Landau, Verwey 
and Overbeek (DLVO) presented a theory which accounts 
for the interaction between weakly charged homogeneous 
surfaces^. The net interaction between two surfaces was 
attributed to the electrostatic double layer forces and the 
van der Waals force. The van der Waals force dominates 
when the separation between the surfaces is small, while 
the electrostatic repulsion is dominant on larger length 
scales. The theory works reasonably well for weakly 
charged homogeneously surfaces^ and has been widely 
used to study colloidal stability. It fails, however, to ac¬ 
count for the correlation induced attraction between like- 
charged objects inside an electrolyte solution containing 
multivalent counterions^^— or for ionic specificity^^—. 

It is also well known that two surfaces with annealed 
positive and negatively charged domains feel attrac- 
tion^i^— . In this case, positive domains on one surface 
become correlated with the negatively charged domains 
on another surface, resulting in an attractive interac- 
tioii2S,i2I. Jho et al^ carried out numerical simulations 
for flat surfaces with movable charged domains. Long- 
range attractive force was observed and the mechanism 
behind the attraction was found to be the positional cor¬ 
relation between oppositely charged domains. 

Recently, Silbert et ali^ conducted an interesting ex¬ 
periment to explore the interaction between heteroge¬ 
neously charged surfaces. Remarkably, they observed an 
attraction which extended up to 500 A. At first the at¬ 
traction was attributed to the correlation between the 
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oppositely charged domains on the two surfaces. How¬ 
ever when a rapid shear motion was introduced between 
the surfaces to frustrate the correlations, the attraction 
persisted. Since the charged domains under the shear 
motion could not become correlated, the distribution of 
surface charge was effectively quenched. It became clear 
that the correlations weren’t the mechanism responsible 
for the attraction in the experiments of Silbert et al. 
The question then became: what was? Silbert et al. at¬ 
tributed the attraction to the unequal nature of the repul¬ 
sive and attractive interactions between like-charged and 
oppositely-charged domains. This, however, appeared to 
contradic the conclusion that quenched charge disorder 
should not lead to attraction between heterogeneously 
charged surfaces^Sri^^. To justify their conclusion, Silbert 
et al. presented a simple argument based on the Poisson- 
Boltzmann (PB) equation. They suggested that the in¬ 
teraction between two neutral surfaces with a quenched 
charge disorder arises from an asymmetry in the interac¬ 
tion between like and oppositely-charged domains inside 
an electrolyte solution. For like-charged domains, the 
counterions are required to stay between the surfaces to 
preserve the local charge neutrality, while for the oppo¬ 
sitely charged domains this is not necessary. The entropic 
contribution to the overall force is, therefore, asymmet¬ 
ric in the two cases. Silbert et al., then, suggested that 
the force between the two heterogeneously charged ran¬ 
dom surfaces can be estimated as an arithmetic average 
of the force between two like-charged and two-oppositely 
charged homogeneous surfaces. 

In the present work, we will use Monte Carlo sim¬ 
ulations to show that the conclusions of Silbert et al. 
are qualitatively correct. The long-range attraction be¬ 
tween heterogeneous surfaces with a quenched disorder 
arises due to the asymmetric interaction between oppo¬ 
sitely and like-charged domains. On the other hand, we 
will demonstrate that a simple arithmetic average of the 
force between like-charged and oppositely-charged homo¬ 
geneous surfaces is not sufficient to quantitatively ac¬ 
count for the range and strength of the attraction be¬ 
tween heterogeneous surfaces with a quenched charge 
disorder, and a more sophisticated calculation must be 
performed. The paper is organized as follows. In sec- 
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tion im we explain the model and the simulation details. 
In section m we summarize our results. In section m 
we conclude our work. 


II. THE MODEL AND SIMULATION DETAILS 

In the experiments of Silbert et al. the charged do¬ 
mains were produced by the adsorption of cationic mi¬ 
celles to an anionic substrate. To simplify the calcula¬ 
tions we will neglect the spatial extent of micelles and 
project all the charge onto a flat surface, see Fig. [TJ As 
an additional simplification, we will divide our surfaces 
into positive and negative domains of the same surface 
area. 



+ + 
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FIG. 1. Representation of the domains. Negative sites with 
absorbed surfactants can be considered as positive sites. In 
“a” a side view and in “b” a top view of the wall. 

Our system then consists of two flat surfaces of di¬ 
mensions Lx and Ly, enclosing an electrolyte solution. 
For simplicity we set Lx = Ly = 400 A. The plates are 
separated by a distance L. The solvent is assumed to 
be a uniform dielectric of permittivity e^,. The Bjerrum 
length, defined as Xb = /ksTcw — where q, ks and T 
are the proton charge, the Boltzmann constant, and the 
temperature, respectively — is 7.2 A for water at room 
temperature. The ions are modeled as hard spheres of ra¬ 
dius 2 A. To perform simulations we use a grand canoni¬ 
cal Monte Carlo algorithm (GCMC)^“— , see appendix lAl 
for details. The system is in contact with a salt reservoir 
at concentration pg. As an input, the GCMC requires 
the chemical potential of the ions of reservoir. For 1:1 
electrolyte this can be calculated from ps using the mean 



FIG. 2. Ionic density profiles between two equaly charged 
plates, with charge density —0.01602 C/m^. Symbols repre¬ 
sent simulations data while lines represent PB curves. The 
concentration of the monovalent salt in reservoir is 20 mM. 
The solid line and circles represent positive ions while the 
dashed line and squares, negative ones, a is the position be¬ 
tween two surfaces. 


spherical approximation (MSA)2S“— , which is very ac¬ 
curate for weakly interacting ions. Similarly, MSA also 
provides us with the osmotic pressure of the bulk elec¬ 
trolyte. The force per unit area between the two surfaces 
is then the pressure between the two plates, minus the 
pressure of the bulk (reservoir) electrolyte. The pressure 
on each plate is calculated by taking into account the 
electrostatic interactions and the entropic force arising 
from the momentum transfer during the collisions of the 
ions with the surfaces. The entropic contribution is cal¬ 
culated using the method of Wu et The details 

of the calculation are presented in appendix |B] 

To obtain the chemical potential and the osmotic pres¬ 
sure of a reservoir containing 2:1 electrolyte we first per¬ 
form a bulk GGMG simulation. In this simulation the 
chemical potential of electrolyte is fixed and the aver¬ 
age concentration of ions inside the reservoir is calcu¬ 
lated. We then perform a NPT simulation to calculate 
the osmotic pressure of the electrolyte^^ at this concen¬ 
tration. The NPT MG simulation method is described in 
the appendix O To calculate the electrostatic energy, a 
3D Ewald summation method with a correction for the 
slab geometry of Yeh and Berkowita^ is used. 

Before considering the interaction between heteroge¬ 
neously charged surfaces, we investigate two simpler 
cases: equally-charged and oppositely-charged homeoge- 
neous surfaces, both in contact with a monovalent salt 
reservoir. In this case each plate is formed by point 
charge pseudo-particles, uniformly distributed, with sep¬ 
aration Lx/Ng along the surface. The charge of the 
pseudo particles is adjusted to obtain the desired sur¬ 
face charge density. The system has periodic boundary 
condition in x and y directions. We set Ng = 40. To 
test the simulations, we compare our results with the 
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FIG. 3. Osmotic pressure for two equaly charged plates. The 
parameters are the same as Fig. O 


solution of the Poisson-Boltzmann (PB) equation. For 
weakly charged homogeneous surfaces inside a dilute 1:1 
electrolyte PB equation is expected to be very accurate^. 
The algorithm to solve the PB is the same as in Ref.— 
adapted to the slab geometry. In Fig. [51 we compare 
the density profiles obtained using the simulations to 
the solution of PB equation. As expected, a very good 
agreement between simulations and theory is obtained. 
For homogeneously charged plates the force between the 
surfaces can be easily calculated using the contact the- 
orem^“— . The net force per unit area is the difference 
between internal and external pressures. As mentioned 
above, the osmotic pressure of the reservoir can be ob¬ 
tained using MSA, which agrees perfectly with the NPT 
simulations. On the other hand, within the PB approx¬ 
imation the electrostatic correlations between the ions 
are completely ignored, and the bulk pressure of 1:1 elec¬ 
trolyte is simply that of an ideal gas. In spite of this very 
crude approximation, for the parameters considered we 
see an excellent agreement between the simulations and 
theory, see Fig. [31 In Figs.|l|and|Sl we show the ionic dis¬ 
tribution and the force per unit area for two oppositely 
charged homogeneous surfaces. Again the agreement be¬ 
tween GCMC simulations and PB equation is very good. 

To calculate the force between two randomly charged 
heterogeneous surfaces we divided each plate into equi- 
sized domains, half of which are positively charged while 
the other half are negatively charged. For each charge 
distribution, we calculate the force between the two sur¬ 
faces. The net force is then calculated as an arithmetic 
average over the quenched disorder. Because the number 
of configurations grows exponentially fast with the num¬ 
ber of charged domains, in this paper we will consider 
surfaces with only two and four regions. 

We start by considering two overall neutral plates with 
two charged domains each: one positive and one neg¬ 
ative, as illustrated in Fig. [51 There are two possible 
configurations, Al-Al and A1-A2. The net force is the 


FIG. 4. Ionic density profiles between two oppositely charged 
plates, with charge density ±0.01602 C/m^. Symbols repre¬ 
sent simulations data while lines represent PB curves. The 
concentration of the monovalent salt in reservoir is 20 mM. 
The solid line and circles represent positive ions while the 
dashed line and squares, negative ones, a is the position be¬ 
tween two surfaces. 



FIG. 5. Osmotic pressure for two oppositely charged plates. 
The parameters are the same as Fig. (4) 


average over these two configuration. We next consider 
the interaction between two plates containing 4 charged 
regions, two positive and two negative. This is illustrated 
in Fig. |7l There are 6 possibilities for each plate. As a 
result, for two interacting plates, there are 36 distinct 
configurations. However, various of these configurations 
are degenerate and are connected by symmetry. We cal¬ 
culated the energy for each configuration and obtained 
the degeneracy factors which are presented in Tab. HI All 
parameters for this model are the same as in the pre¬ 
vious case. By dividing the plate into four areas, the 
positive and negative domains become smaller, in com¬ 
parison with the two-region model. This way we are able 
to explore the effect of domain size on the overall in¬ 
teraction between the surfaces. Furthermore, note that 
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FIG. 6. Representation of surfaces when they are divided into 
two charged domains. The surfaces are overall neutral. 


because of the periodic boundary conditions imposed by 
the Ewald summation, the plates are actually infinite, 
with a charge distribution corresponding to stripes or a 
checkerboard patterns. 


Configurations: Bl-Bl B1-B2 B1-B3 B1-B5 B1-B4 B3-B3 
Degeneracies: 2 2 16 8 4 4 


TABLE I. Number of different configurations for the four- 
region model. 
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FIG. 7. Representation of surfaces when they are divided into 
fonr charged domains. The snrfaces are overall nentral. 


III. DISCUSSION 
A. Monovalent salt ions 



FIG. 8. Net force per unit area between two neutral sur¬ 
faces with absolute site charge density equal to 0.056 C/rr?, 
in contact with a salt reservoir at concentration 0.01 M. The 
circles and squares represent two and four-region models, re¬ 
spectively. The lines are guides to the eye. 


and the entropic force is reduced. Calculating the aver¬ 
age over the two charge distributions, we see that the 
net force is repulsive at short distances, but becomes at¬ 
tractive at larger separations, see in Fig. [S] It should 
be mentioned that in all cases the electrostatic contribu¬ 
tion to the force is attractive, as is expected for a charge 
neutral system. The net repulsion at short distances is 
observed because of a strong entropic force produced by 
the confined counterions. 

In Fig. [51 we plot the net force between two randomly 
charged surfaces for various surface charge densities of 
the two domains. As can be seen, by reducing the charge 
density, the attraction between the plates decreases. This 
is not surprising since the attraction is caused by the 
electrostatic interaction. 

To investigate the effect of domain size on the inter¬ 
action between the surfaces, we divided the area of each 
plate into four regions. By doing this, each charged do¬ 
main becomes smaller. There are a total of 36 different 
configurations many of which, however, are related by 
symmetry. We find that there are only 6 distinct ar¬ 
rangements, each with its own degeneracy factor listed 
in Tab. H The net force can now be easily calculated as a 
weighed average over the configurations listed in Tab. [H 
Similar to the two-region model, a net attraction is ob¬ 
served, see Fig. [51 We see, however, that the attraction 
is somewhat weaker than for the two-region model. 


In the two-region model, the plates are divided into 
two domains, one positive and one negative, as shown in 
Fig. m The surface charge densities of both domains are 
the same. In Fig. [51 we plot the average (over disorder) 
force per unit area between the two plates. For con¬ 
figuration Al-Al, at all distances repulsion is observed. 
In configuration AI-A2, the two plates feel attraction. In 
this case the number of ions between the plates is smaller. 


B. Divalent salt ions 

We next explore the effect of the charge asymmetry 
of electrolyte on the interaction between two heteroge¬ 
neously charged surfaces. We consider a reservoir of 2:1 
electrolyte at concentration 10 mM. The surface patterns 
and surface charge densities are the same as for the mono- 
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FIG. 9. Net force per unit area between surfaces with 
two charged region. The reservoir 1:1 salt concentration is 
0.01 M. The circles and squares represent the charge densi¬ 
ties 0.04 C/re? and 0.072 C/m^, respectively. The lines are 
guides to the eye. 


valent salt. Fig. [TU] shows that a charge asymmetric elec¬ 
trolyte leads to a stronger attraction between two het¬ 
erogeneous random surfaces than a 1:1 electrolyte. Fur¬ 
thermore, in the case of 2:1 salt, the difference between 
the force for two and four-region models is larger, see 
Fig. [ini showing that the size of the charged domains is 
more important for the interaction between the surfaces 
in asymmetric electrolyte solutions. In Fig. [TTl we plot 
the net force for various surface charge densities. Simi¬ 
lar to what was observed for 1:1 electrolyte, decreasing 
the surface charge density of the domains diminishes the 
attraction between the plates. 



FIG. 10. A net force per unit area between two neutral sur¬ 
faces with absolute site charge density equal to 0.056 C/m^, 
for 2:1 salt reservoir at concentration 10 mM. The circles and 
squares represent two and four sites model, respectively. The 
lines are guides to the eye. 



FIG. 11. Net force for various surface charge densities of the 
domains. The two-region model is used. The reservoir 2:1 salt 
concentration is 0.01 M. The circles and squares represent the 
charge densities 0.04 C/rr? and 0.072 C/m^, respectively. The 
lines are guides to the eye. 



FIG. 12. Comparison of the calculation of Silber et ah— with 
the results of the GGMC simulation for the two-region model 
for various charge densities of the domains, from top to bot¬ 
tom 0.04, 0,056 and 0.072C'/m^. The calculation of Silber 
et al. significantly overestimates the attraction. The devia¬ 
tions grow with the increasing surface charge density of the 
patches. Solid lines are the results of GGMC while the dashed 
lines are the results of the model of Silber et al. 


IV. CONCLUSION 

We have presented a simple model of interaction be¬ 
tween randomly charged heterogeneous surfaces inside an 
electrolyte solution. Two and four-region models were 
explored. Unfortunately due to the exponential growth 
of configurations, a study of surfaces with smaller charge 
domains is not viable. For example a surface with 16 do¬ 
mains would lead to a total of 165636900 different con¬ 
figurations, a brute force study of which is clearly im¬ 
possible. Nevertheless, exploration of two and 4 region 
models has already provided a number of valuable in- 
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sights. Indeed as was argued by Silbert et al., interac¬ 
tion between two overall neutral surfaces with randomly 
charged domains is attractive at large distances. We find 
that the attraction decreases with the size of domains and 
increases with the charge asymmetry of the electrolyte 
solution. Silber et al. attributed the attraction between 
neutral randomly charged surfaces to the asymmetry of 
the interaction between like and oppositely-charged do¬ 
mains. While the interaction between like-charged do¬ 
mains has a strong entropic component, the entropy is 
less important for oppositely charged domains, since the 
counterions are not required to stay between the oppo¬ 
sitely charged regions to neutralize their charge. Based 
on this observation, Silber et al. concluded that the inter¬ 
action between two heterogeneous random surfaces can 
be estimated as an arithmetic average of the force be¬ 
tween two infinitely large like-charged surfaces and the 
force between two oppositely-charged surfaces. As we 
saw in Section m the interaction between like-charged 
and oppositely charged surfaces can be very accurately 
calculated using the PB theory. We can, therefore, easily 
check the model of Silber et al. by comparing it with our 
simulations. In Fig. [T^] we contrast the model of Silber 
et al. with the results of our GCMC simulations for the 
two-region model. Although qualitatively correct, the 
calculation of Silber et al. significantly overestimates the 
attraction between the two surfaces. The error increases 
with the increasing surface charge density of the charged 
domains. Finally, we expect that in the limit that the 
area of domains becomes small, the approach of Naji and 
Podgornik for randomly charged surfaces should become 
accurate^. Unfortunately, we are not able to check this 
limit in our model, since, as mentioned previously, reduc¬ 
tion of the domain size leads to an exponential growth of 
configurations. 

Up to now we have only considered overall neutral sur¬ 
faces. In the future it will be interesting two explore the 
effect of breaking the charge neutrality on the attraction 
between heterogeneous randomly charged objects. 
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Appendix A: Grand Canonical Monte Carlo Simulation for 
a:l salt 


state i is proportional to^^ 


YiN++N.)^-l3Ei+l3N+tJ.+ +l3N.^i. 

N+\N_\A^_^+ A^_^- 


(Al) 


where V is the volume accessible to particles, N± are 
the number of cations and anions, Ei is the electrostatic 
energy of the state, /i± are the chemical potentials and 
A± are the thermal de Broglie wavelengths. 

The transition probability for addition is: 


p. ya+l^-PEj+PEi+PiJ.^+aPii- 

^ {N+ + 1)(A^_ + a){N_ + a - 1)...(A^_ + l)A^Ai“ ’ 

(A2) 

where j is the state after addition. We define the parame¬ 
ter z = /A^A^“. This parameter is the input 

of simulations. The transition probability for addition 
can be rewritten as 

p. ^ z ya+l^-PE^+pEi 

^ “ {N+ + 1){N_ + a){N_ +a- 1)...{N_ + 1) ’ 

(A3) 

The thansition probability for the removal is: 

p, - l)...(jV_ - g + 1) 

p, z U“+i ’ ^ ’ 

where j is the state after removal. After the transition 
probabilities have been calculated, they are compared 
with a random number, uniformly distributed between 
0 and 1. If this random number is lower than the transi¬ 
tion probability, the movement is accepted. Otherwise it 
is rejected. 


Appendix B: Entropic force 

The algorithm of Wu et alM, constructed to calculate 
the entropic force between two colloidal particles, can be 
easily adapted to the planar geometry. The entropic force 
is give by the expression 


^F = 


(Ni) 

Az 


(Bl) 


where A^i is the number of overlaps of one of the walls 
with the free ions (which are held fixed) after a displace¬ 
ment Az. The force is obtained in the limit Az 0. 
The entropic pressure is 


PP 


(Ni) 

AzLxLy 


(B2) 


Here, we briefly review the GCMC method for a:l salt. 
In each step, we have three possibilities, simple move¬ 
ment of ions and the addition or removal of ions. In 
order to keep the charge neutrality, if one randomly adds 
or removes a cation with valency a, a anions must also 
be added or removed. The probability of a particular 


Appendix C: NPT Monte Carlo simulations 

In this appendix, the NPT MC simulation is rapidly 
reviewed. Besides the particle movement, the volume of 
the simulation box is varied. This kind of movement is 
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randomly chosen with a small probability, normally, 1 /N, 
where N is the number of particles in the box. After the 
increment or decrement of the volume of the box, the 
particle positions are rescaled. The transition probability 
of acceptance of the changes in volume is 

Pi = (Yl\ ^-PEj+pEi-PiVj-Vi) ^ 

Pi ) 

where j is the new state, Vi and Vj are the volume ac¬ 
cessible to the N particles in the old and new states, 
respectively. The pressure P is the input of the simula¬ 
tion. Again, if a random number uniformly distributed 
between 0 and 1 is lower than the transition probability, 
the movement is accepted. Otherwise it is rejected. 
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